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Abstract. A Poisson-Nernst-Planck-Fermi (PNPF) theory is developed for studying 
ionic transport through biological ion channels. Our goal is to deal with the hnite size of 
particle using a Fermi like distribution without calculating the forces between the particles, 
because they are both expensive and tricky to compute. We include the steric effect of ions 
and water molecules with nonuniform sizes and interstitial voids, the correlation effect of 
crowded ions with different valences, and the screening effect of water molecules in an inho¬ 
mogeneous aqueous electrolyte. Including the hnite volume of water and the voids between 
particles is an important new part of the theory presented here. Fermi like distributions 
of all particle species are derived from the volume exclusion of classical particles. Volume 
exclusion and the resulting saturation phenomena are especially important to describe the 
binding and permeation mechanisms of ions in a narrow channel pore. The Gibbs free energy 
of the Fermi distribution reduces to that of a Boltzmann distribution when these effects are 
not considered. The classical Gibbs entropy is extended to a new entropy form — called 
Gibbs-Fermi entropy — that describes mixing conhgurations of all hnite size particles and 
voids in a thermodynamic system where microstates do not have equal probabilities. The 
PNPF model describes the dynamic how of ions, water molecules, as well as voids with 
electric helds and protein charges. The model also provides a quantitative mean-held de¬ 
scription of the charge/space competition mechanism of particles within the highly charged 
and crowded channel pore. The PNPF results are in good accord with experimental currents 
recorded in a 10®-fold range of Ga^"*" concentrations. The results illustrate the anomalous 
mole fraction ehect, a signature of L-type calcium channels. Moreover, numerical results 
concerning water density, dielectric permittivity, void volume, and steric energy provide use¬ 
ful details to study a variety of physical mechanisms ranging from binding, to permeation, 
blocking, hexibility, and charge/space competition of the channel. 


1. INTRODUCTION 

Biological functions of proteins depend on the details of the mixtures of ionic solutions 
found outside and inside cells. Trace concentrations (< 10“® M) of calcium ions (Ga^^) 
and other signaling molecules provide physiological control of many biological pathways and 
proteins inside cells [ij. For example, voltage-gated calcium (Gay) channels exhibit the 
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anomalous mole fraction effect that effectively blocks abundant monovalent cations by a 
trace concentration of Ca^’*' ions The fundamental mechanism of the calcium channel 

is of great technological and biological interest 0,0. Multiscale analysis seems necessary 
since calibrated all atom simulations of trace concentrations of ions in physiological solutions 
are not likely to be available in the near future. 

Interactions between diffusion and migration in the electric held are central to the biol¬ 
ogists’ view of channels 0,0. Following the drift-diffusion (DD) model in semiconductors, 
Eisenberg et ah have proposed the Poisson-Nernst-Planck (PNP) model to calcu¬ 

late rather than assume 0, 0 the electric held and then the ionic current in biological ion 
channels. Interactions of ions and hows in narrow channels, and saturation in binding sites 
are also central to the biologists’ view of channels 0, 0. It has been difficult to combine 
the two views — computed electric helds in channels showing interactions and saturation 
— because charges in PNP (and quasiparticles in DD) are points with no diameter and so 
cannot saturate the aqueous channel [ij] the way real ions do. Recently, the saturation of 


spheres of any size has been described by a Fermi like distribution derived in [12| from the 


conhguration entropy of mixtures of ions of any diameter and composition. The steric ehect 
has been shown to be very important to adequately describe equilibrium systems [iMi. 

We extend the Poisson-Fermi model [l0 in two important rather novel ways. We include 
the excluded volume of water molecules and the ‘empty space’ created by packing constraints 
and voids between particles. The equilibrium model is also generalized here to a nonequilib- 
rium model called the Poisson-Nernst-Planck-Fermi (PNPF) model that can describe flow, 
including the steric effect of all particles, the correlation effect of ions and water molecules, 
the screening effect of water, as well as the charge/space mechanism in the channel pore 
at and away from equilibrium. This treatment unites diffusion and electric current, with 
interactions and binding in narrow channels. 

Both discrete and continuum forms of Gibbs free energy of the electrolyte are developed 
in this paper. The Gibbs-Fermi free energy functional for the Fermi distribution is shown 
to reduce to the Gibbs-Boltzmann functional for the classical Boltzmann distribution when 
both steric and correlation effects are not present. Moreover, a new entropy form called the 
Gibbs-Fermi entropy is proposed here to connect the spatial distribution of ions, water, and 
voids between them (that may vary) with the change of local probabilities of each species 
(and void volume) which of course usually have different sizes. The Gibbs-Fermi entropy 
is a consistent generalization of the global Boltzmann entropy and the classical local Gibbs 
entropy widely used to describe systems without steric and packing constraints. 

The steric effect. The steric effect of crowding produces a steric energy term in PNPF 
that is a quantitative statement of the crowded charge effect of charge/space competition. 
The charge/space competition theory introduced by Nonner and Eisenberg to explain cal¬ 
cium selectivity has been developed in a long series of papers using Monte Garlo methods 
by Boda and Henderson, and density functional methods by Gillespie and collaborators. 
This ‘all spheres’ approach successfully describes almost all selectivity properties of calcium 
channels and the main properties of sodium channels such as the micromolar Ga^^ affinity 
for L-type calcium channels [l5-T3|, the wide range of Ga^^ affinities for different types of 
calcium channels, and the switch in selectivity from calcium to sodium when the side chains 
of the selectivity hlter are switched from EEEE (gin gin gin gin) to DEKA (asp gin lys 
ala) IlSl. Il8|. It also accounts for the selectivity between monovalent cations of different size 
19-2l|| and for the self-organized pore structures for selective ions [10. In the biologically 
crucial and special case of Na^ vs K"*" selectivity in the DEKA sodium channel (so central 
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to the function of the nervous system and metabolic budget of mammals with large brains 
(^. 2^), control variables can even be identihed that independently control selectivity and 


binding [l^, Il8 


Interactions. Our main goal is to show how interactions of diffusion, electrophoretic 
migration, steric exclusion, and imperfect packing of particles can all be treated quantita¬ 
tively in a unihed framework to analyze binding and flow in crowded ion channels, without 
explicitly calculating forces between individual ions, water, or voids. We show that a Fermi 
like distribution is able to describe these interactions well enough to account for a wide range 
of important properties of ions in channels. We wonder how well this approach can describe 
the myriad of nonideal properties reported in the physical chemistry experimental literature 
(for more than a century) which have escaped canonical description up to now 

Numerical results produced by the PNPF model are in accord with the experimental 
data reported by Aimers and McCleskey in 1984 for Cay channels over a 10®-fold range of 
concentrations of calcium ions [^. Their experimental data has been a benchmark for selec¬ 
tivity ever since. Their data has been used as a tareet for models using a variety of methods 
ranging from physiological and crystallographic to molecular dynamics (MD) 


[3a, IM 


Monte Carlo (MC) [ISL Il7l. |40| . as well as continuum 


Brownian dynamics (BD) [3^39 
approaches [411^4^. 

The remaining part of the paper is organized as follows. A derivation of the conhguration 
entropy of all hard-sphere ions and water with voids is proposed in Section 2, where a 
Fermi type of excess chemical potential, Gibbs-Fermi free energy functional, and Gibbs- 
Fermi entropy are also introduced. All these models seem to be new to the literature, 
as far as we know, because they treat hnite size water molecules and voids explicitly. In 
Section 3, we extend the Gibbs-Fermi theory to the Poisson-Nernst-Planck-Fermi theory for 
studying ionic transport, steric energy, water density, and void distribution in equilibrium 
or nonequilibrium conditions. In Section 4, a molecular-continuum model specihc to L-type 
calcium channels is presented to show how to implement in a consistent way the PNPF 
theory of the molecular hlter region of few particles joined to the bath region of numerous 
particles. Section 5 demonstrates that PNPF currents agree quite well with the experimental 
currents reported in under the same membrane potential and the same 10®-fold range 
of Ga^^ concentrations measured in the experiment. These conductance results seem to £t 
data better than other models we know of. Some concluding remarks are given in Section 6. 


2. FERMI DISTRIBUTION AND GIBBS-FERMI ENTROPY 


Based on the configurational entropy model proposed in [12[ for aqueous electrolytes with 
arbitrary K species of nonuniform size, hard spherical ions, we extend the free energy of the 
model to 


K+l 


F{N) = (j)J2 ^ 


( 1 ) 




by including specifically the excluded volume effect of the next species [K + 1) of water 
molecules. Here, 0 is the electrostatic potential, Nj is the total number of j species particles 
carrying the charge Qj = ZjC with the valence Zj, e is the proton charge, is the Boltzmann 
constant, and T is the absolute temperature. The volume of a j type particle is Vj = 
with radius aj. It is important to note that water is treated as a polarizable hard sphere 
with zero net charge in Eq. (1), so zk+i = Qk+i = 0. The polarizability of water and the 
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inclusion of voids represent important generalizations from the classical primitive solvent 
model used to describe calcium channels ji^. The last term in (1) describes the mixing 
entropy of all ions and water molecules over a total of N available nonuniform sites in a 
system with 

^ N\ 

= n ^ ■ (2) 


i=i 


(nfjw^!) (A'-Ef.ViV,)! 


where Wi = N\/{Ni\{N — Ni)\) is the number of combinations for the distribution of Ni in 
all vacant sites N. W 2 = {N — Ni)\/{N 2 l{N — A^i — N 2 )\) is the number of combinations 
for the distribution of N 2 in N — Ni vacant sites after Ni being distributed, and so on. 
After all particles are distributed, there remains (in this model) just a single site Nk +2 = 
N — ^1,^=1 1 used to represent the (continuously connected) voids created by 

defects in the packing structure of all particles of all types and by Coulomb and steric forces 
(e.g., Lennard-Jones) between particles. This void structure is represented as the last species 
K + 2 in our model. We are unaware of other all-spheres models that deal explicitly with 
the voids between spheres. We suspect that including such voids is needed because voids 
are in different amounts depending on the composition of the solution and can move in any 
system of spheres crowded into a small space. 

Obviously, all properties of water cannot be represented this way: water is a highly 
charged molecule although its net charge is zero, and polymeric structures can exist and 
may be important, along with hydrogen bonds of the low or high energy type (dsl. 46 


Moreover, not all defects in packing can be represented by a single void site, just as not 
all properties of water can be represented by uncharged spheres. The question is whether a 
model that includes only the excluded volume of water and a continuous void space between 
particles is able to deal with the selectivity data of the calcium channel. We will see that it 


can. 


The total volume V of the system consists of the volumes of all particles and the total 

Under the bulk condition, dividing this 


void volume vx+ 2 , he., V = Yl!j=i 
equation by V yields the bulk void volume fraction 


VjNj Vk+2- 


r^ = 


'Vk+2 

V 


= 1 


u 


i=i 


= 1 


K+l 


z. '’iU 

i=i 


(3) 


expressed in terms of the nonuniform particle volumes Vj and the bulk concentrations (7® 
of all particle species. We are aware that a model of this sort can be extended into an 
all-spheres model of ionic solutions of the so called bio-ions Na"*", K"*", Ca^^, and Cl“. 

Using the Stirling formula InM! M\nM — M with M » 1, the electrochemical 
potential of particle species i = l,---,ih-|-lis 




dF{N) 

ON, 


= qiif) + ksT In 




N - V 


(4) 


from which we deduce global probabilities Pi = Ni/N for all particle species. If we extend 
our theory by introducing local probabilities Pi{r) = UjUj(r) that depend on location, in 
effect allowing probabilities to depend on location as in the theory of stochastic processes 
(applied for example to ionic channels jlll, 49|), the electrochemical potential can 
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be generalized locally to 

/ii(r) = (?i0(r) + fcsTln 

v^Cf 


ViCi{r) 


= gi0(r) + fcsTln + //“(r) 


/if(r) = ksTln 


r(r) 


K+l 

’ r(r) = 1 - ^ ^'iC'J(r) = Vk+ 2 Ck+ 2 {j:), 


(5) 

( 6 ) 


i=i 


where C'i(r) is the concentration function of spatial variable r in the solvent domain 
/i®^(r) is the excess chemical potential, and r(r) is the void fraction function with Ck+ 2 {'^) 
representing the distribution function of interstitial voids. When 0 = 0, C'i(r) = Cf and 
hence = /rf = fc^Tln (ujCf/r®) is a constant. 

The excess chemical potential is a measure of nonideality that helps understand qual¬ 
itative behavior. For example, the larger the size Vi of a type i particle, the larger is the 
activation barrier and the harder it is for the particle to make a transition at r from 

a local minimum of /i|^ to another local minimum nearby [s^. The transition mechanism 
is related to the vacancy conhguration as well, i.e., the smaller the value of r(r), the more 
crowded the ions are at r, the harder transition. The excess chemical potential is closely 
related to the sizes of all particles their interstitial voids r(r), their conhgurations C'j(r), 
as well as their bulk concentrations . 

To our knowledge, all existing continuum models do not explicitly take into account the 
hnite size effect of water, let alone the effect of interstitial voids. We did not consider these 
two effects in our previous work 1^ in which water was treated as a single continuous 
dielectric medium without any voids and the resulting electrochemical potential /ij(r) was 
shown to be a mathematical description of the primitive model of electrolytes, as used in 
most Monte Carlo and density functional theory models. Our continuum primitive model 
could well match Monte Carlo (discrete primitive model) results that were obtained in 
equilibrium state. However, as we proceeded to study nonequilibrium systems using this 
primitive model, we had difficulty computing the experimental currents reported in [2| due 
to either inconsistent physics or divergent numerics. 

The calcium channel operates very delicately in physiological and experimental condi¬ 
tions as it shifts its exquisitely tuned conductance from Na’''-flow, to Na+-blockage, and to 
Ca^'''-flow when bath Ca^"*" concentration varies from 10“^° to 10“^ M. The 10®-fold range 
of experimental conditions make modeling and numerical implementation very challenging. 
This huge dynamic range was accommodated in our previous work by using an artihcial 
potential to conhne mobile oxygen ions of side chains within a hlter region, just as that used 
in all Monte Carlo simulations on the same channel (see e.g. 18|). The artihcial potential 


hindered our effort to match the experimental data since it is a gross approximation of the 
constraining energy needed to keep the protein atoms in the hlter region without specihcally 
considering the void ehect. Indeed, using a restraining potential can lead to inconsistencies, 
since maintaining the steric and electrical potential as conditions change requires injection 


of energy and charge into the system . We obtain convergent and consistent results using 
the steric potential in place of the artihcial constraining potential of earlier models. The 
steric potential is an output of our model and varies automatically as conditions change. It 
has the same units as a conhning potential but is as diherent as the voltages at an input 
and an output of an ideal ampliher. The following analysis shows that the void species in 
our model is important not only to describe a consistent physics of the steric potential but 
also to compute the steric energy that can rehect the 10®-fold experimental conditions. It 
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will be interesting to examine the properties of a model of bulk ionic solutions that contains 
voids calculated consistently in an all-spheres model of ions and water. 

Setting /ii(r) = /if (see below for physical reason), the concentration of species i particles 
can be expressed by the Fermi like distribution function 

C,(r) = Cfexp(-/?,.#.(r) + S“(r)). S“(r) = liih^, (7) 

where (3i = qi/ksT and is called the steric potential that describes the combined effect 

of all excess chemical potentials /i®^ of all particle species j = 1, • • ■ iF -|- 1. The distribution 
(7) is of Fermi type since all concentration functions 


CAv) = 


Oti 


E r 
i 




1 -h aiVi 


< 


(l/ai) + Vi 


1 

Vi 


( 8 ) 


i = 1, • • • , iF -|- 1, are bounded from above with a* = Cf exp (—/9j0(r)) /F® > 0, i.e., C'i(r) 
cannot exceed the maximum value 1/vi for any arbitrary (or even inhnite) potential 0(r) at 
any location r in the domain In this mean-held Fermi distribution, it is impossible for 
a volume Vi to be completely hlled with particles, i.e., it is impossible to have ViCAv) = 1 
(and thus F(r) = 0) since that would make the excess chemical potential /i|^ inhnitely large 
or = —oo and hence C'i(r) = 0, a contradiction. For this reason, we must include 
the voids as a separate species if water is treated as hard spheres. Otherwise, 
the volume Vi would be easily hlled by particles in the mean-held sense at moderate electric 
potential such that the steric potential would be unphysical. The requirement of voids when 
all particles are represented as hard spheres will be justihed again from a viewpoint of Gibbs’ 
free energy. 

The classical Boltzmann distribution appears if all particles are treated as volumeless 
points, i.e., Uj = 0 and F(r) = F® = 1. It may produce an inhnite concentration Gj(r) —)■ oo 
in crowded conditions when —/5j0(r) —)■ oo, close to charged surfaces for example, an im¬ 
possible result 12-1^. The difficulty in the application of classical Boltzmann distributions 
to saturating systems has been avoided in the physiological literature (apparently starting 
with Hodgkin, Huxley, and Katz by redehning the Boltzmann distribution to deal with 
systems that can only exist in two states. This redehnition has been vital to physiological 
research and is used in hundreds of papers [E^, but confusion results when the physiolo¬ 
gists’ saturating two-state Boltzma nn is not kept distinct from the unsaturating Boltzmann 
distribution of statistical mechanics 


55 


To further account for the correlation ehect of ions and the screening ehect of water 


molecules, we have developed efficient 3D methods [13| for solving the Poisson-Fermi (PF) 
equation [i2-[Ii,H[53 


K 


- 1) VV(r) = ?.Ci(r) = p(r) 


(9) 


i=l 


self-consistently with Eq. (7) for 0(r), where A is a correlation length [5g, [STj, = e„eo, 
Cw is a dielectric constant of water in the bath, and Cq is the vacuum permittivity. The 
fourth-order PF equation reduces to the classical Poisson-Boltzmann (PB) equation when 
Ic = S'*’'®(r) = 0. If Ic 7^ 0, the dielectric operator F = 6^(1 — is used to approximate the 
permittivity of the bulk solvent and the linear response of correlated ions js^ . The dielectric 


















7 


function ?(r) = es/(l + rj/p) is a further approximation of e". It is found by transforming 
Eq. (9) into two second-order equations — 1) d' = p and = 'h. We introduce 

a density like variable d' that yields a polarization charge density rj = — p using 

Maxwell’s hrst equation 0(3. 

The free energy formula (1) is useful for a thermodynamic system that involves a limited 
number of particles for MD or MC simulations particularly without flow, or spatially nonuni¬ 
form boundary conditions. If the system is nonequilibrium or has numerous particles and 
complicated boundary conditions, the PF equation (9) will be more suitable for theoretical 
investigation. 

Free energy functional. We look at our model now from the perspective of a gener¬ 
alization of free energy that we call the Gibbs-Fermi free energy. The PF equation is a 
minimizer of the following Gibbs free energy functional 


^Fermi^ / dr 
JUs 




9 = 


/K+l 

PbT ^ 

\i=i L 


Cj In {vjCj) — Cj — Cj In {VK+ 2 CK+ 2 ) + 


ksT 


( 10 ) 


by taking energy variations with respect to cj), i.e., 
follows from = 0, where 


^gFermi 


= 0. The Fermi distribution (7) 


A, = -p? = -iBrin^ (11) 

is the Lagrange multiplier for the mass conservation (the total number W) of particle species 
i jssf- The minimization = 0 is equivalent to setting Pi{r) = pf for Eq. (7) with 

the identity P = vk+ 2 Ck+ 2 - The electrochemical potential (5) can also be dehned by the 
functional as 

h* — -^ hi • (12) 

When If. = S''’'''’’(r) = Vj = 0, j = 1, ■ ■ ■ ,K+1 (without correlation and steric terms), this 
functional yields the PB equation —esV^0(r) = p(r) and the Boltzmann distribution Gj(r) = 
Gf exp (—/dj0(r)) since vk+ 2 Ck +2 = P® = 1 and g = [^j ~ ^j\ ■ Note 

that all Vj in g are canceled before setting Vj = 0 since 

Ci HvA) + ^1)^ = C'j '“('’jFl) - Cj In ^ = C, In (13) 

We need Vj in g to justify that the local electrochemical potential (5) can be dehned by the 
functional (10) and that the Fermi distribution (7) is a consequence of mass conservation 
by (12). 

Voids are needed. To establish a consistent generalization from Boltzmann to Fermi 
distribution it is critical to express the energy functional (10) by means of the void fraction 
function VK+ 2 CK+ 2 {r) = 1 - VjCj{r) in g. Otherwise, In (ui^+ 2 Gi^+ 2 (r)) = Inl = 0 

(without the void term) implies that ~ '^j ~ 0 ^'^'1 Boltzmann 

functional of volumeless particles that we are seeking to replace. This means that it is 

















impossible to treat all ions and water molecnles as hard spheres and at the same 
time achieve a zero volnme of interstitial voids between all particles. Therefore, the 
Gibbs-Fermi functional (10) is not only consistent but also needed (and of course physical) 
with either Vj = 0 (all particles are volumeless points) or Vj ^ 0 (all particles are spheres). 

Of course, we could treat ions as spheres but water as a continuous medium (without 
voids) that then forms the single site in Eq. (2) in place of the voids as previously proposed 
in our paper for the primitive solvent model. The void fraction r(r) would then become 
the water fraction r(r) = 1 —where the upper limit is K instead of iF+1. This 
is precisely the primitive model implemented in the Monte Carlo simulations of Boda and 
Henderson. Consequently, the primitive model may yield incorrect water densities, pressures, 
and dielectric coefficient in mean-field sense for nonideal and inhomogeneous systems (see 
Section 5 for more details). This important limitation in the continuous water version of 
the all-spheres model was pointed out early in its history ji^ ]. 

Comparison with other treatments of finite sized particles. All existing free en- 

include either uniform size (n = Vj for all j) |^, I59l-l62 


ergy functionals that specificall ^ ^ ^ __ 

or nonuniform sizes (n* 7 ^ Vj) |l3|, l 6 ^ cannot reduce to their corresponding Boltzmann 


functionals by directly setting Vj = 0 [60, | 6 ^ because those functionals retain the local 
probability form of VjCj In {vjCj) = pj Inpj in their Gibbs entropy. They use an inconsistent 
reciprocal term involving a uniform particle size, namely 1 /n, instead of a consistent term 
involving the nonuniform particle sizes. The local probability Pj{r) of any particle species j 
in our Gibbs entropy 


K+l 


Cj In (vjCj) - Cj In {VK+ 2 CK+ 2 ) + 
j=i 


keT 


K+l 

-kB Y. a 

i=i 


1 C'.TB 


(14) 


is instead expressed in terms of the local probability ratio Cj(r)/r(r) and the global prob¬ 
ability ratio T^/C® between the particle fraction (probability) C'j(r) and the void fraction 
r(r) per unit volume. In other words, the local probability Pj(r) in the Gibbs-Fermi treat¬ 
ment changes with varying configurations of all particles (F(r) = 1 — 
voids (F(r) = vk+ 2 Ck+ 2 {'^))- The local probability at any location, including the binding 
site, is also connected to the bulk conditions in the bath as implied by Eq. (7). It also 
depends implicitly on the sizes of all particles, valences of ionic species, and long range 
(Goulomb) as well as short range (Lennard-Jones) distances between all particles. All these 
physical properties are lumped into the steric potential functional S'*’''^(r) = In 
simple and unified way. 

The void fractions F® and F represent the Lennard-Jones distances between all particles in 
a mean-field approximation. More specifically, the L-J potential l/(r) = 4 ((cx/r)^^ — (a/r)®) 


m a very 


65l | can be used to determine the distance r between any pair of particles. In bulk solutions. 


the distance r = a yields V{r) =0 that corresponds to a finite but fixed distance a be¬ 
tween any two adjacent particles in the system and thus to the constant bulk void fraction 
F^. Similarly, the nonuniform void function F(r) corresponds to nonuniform inter-particle 
distances r that may or may not equal a for all pairs of adjacent particles. Nonzero L-J 
potentials are in general highly oscillatory and extremely expensive to compute in a system 
of numerous particles. Including external fields adds problems of consistency with spatially 
nonuniform far field boundary conditions to the problems of computational expense. The 
void function F(r) or equivalently the steric potential S'’^’’'^(r) is on the other hand quite 
smooth and relatively very easy to compute. The convolutional density functional on any 
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pair of concentration functions C'i(r) and Cj{r') with a L-J kernel or DFT representation of 


the interaction potential proposed by Eisenberg et ah [6^ is another mean-held approxima¬ 


tion, which is more accurate (but much more difficult to compute reliably) than the steric 
potential since the convolutional functional is nonlocal whereas is local. The 

local steric potential of this paper can be used in place of the nonlocal L-J or DFT potential 
of and therefore the energy variational theory based on the Onsager dissipation principle 
developed in can be applied to the Gibbs-Fermi functional (10). 


3. POISSON-NERNST-PLANCK-FERMI THEORY 


For nonequilibrium systems, the classical PNP model can then be generalized to 

the Poisson-Nernst-Planck-Fermi model by coupling the hux density equation 


dCi{r,t) 

dt 


-V ■ Ji(r,t), r e Ds 


(15) 


of each particle species i = 1, • • • , iP -|- 1 (including water) to the PF equation (9), where 
the hux density is dehned as 


Ji(r, t) = -Di [VGi(r, t) + ^(^^(r, t) V0(r, t) - Ci{r, t) VS'*''‘'(r, t)] , (16) 

Di is the dihusion coefficient, and the time variable t is added to describe the dynamics of 
electric 0(r, t) and steric potentials. The hux equation (15) is called the Nernst- 

Planck-Fermi equation because the Fermi steric potential S'*'''^(r, t) is introduced to the 
classical NP equation. 

At equilibrium, the net how of each particle species is a zero vector, i.e., Ji(r) = 0 (in 
a steady state) which implies that C'jexp(/3j0 — = const. = Cf for 0 = = 0. 

Therefore, the NPF equation (15) reduces to the Fermi distribution (7) at equilibrium. 
Similarly, the classical NP equation reduces to the Boltzmann distribution at equilibrium. 

The steric force. The gradient of the steric potential in (16) represents an en- 

tropic force of vacancies exerted on particles. The negative sign in —means that 
the steric force is in the opposite direction to the ‘dihusional’ force VC*. 

The larger S'*'‘^(r) = In (meaning more space available to the particle as implied by 
the numerator) at r in comparison with that of neighboring locations, the more the entropic 
force pushes the particle to the location r, which is simply the opposite mechanism of the 
diffusional force VGi(r) that pushes the particle away from r if the concentration is larger at 
r than that of neighboring locations. Moreover, the Nernst-Einstein relationship implies 
that the steric flux is greater if the particle is more mobile. Therefore, the gradi¬ 

ents of electric and steric potentials (V0 and VS'*’’^) describe the charge/space competition 
mechanism of particles in a crowded region within a mean-held framework. Since S'*''^(r, f) 
describes the dynamics of void movements, the dynamic crowdedness (pressure) of the how 
system can also be quantihed. 

The motion of water molecules is directly controlled by the steric potential in our model 
and their distributions are expressed by Gi^+i(r, t) = exp (S'*’''^(r, t)). Nevertheless, this 
motion is still implicitly ahected by the electric potential 0(r, t) via the correlated motion of 
ions described by other Cj{r,t) in the void fraction function F(r,t) and hence in the charge 
density p in (9). 
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Water is polarizable in our model. From (9), = \t', and rj = —esd' — p, we 


deduce that the Poisson equation 


e. ^ ^ + r/, 


(17) 


which describes the electric field E = — V0 in the system, contains the charge source not 
only from the ions (p = QiCi) but also from the polar water (p) provided that the 

correlation length Ic ~ IsQi is not zero. The polarization charge density p is proportional to 
the fourth order of the ionic valence Zt. The fourth order dependence shows that even in a 
mean field theory valencw of ions is expected to play an important role as known in chemical 
and biological systems 

In summary, the PNPF model accounts for the steric effect of ions and water molecules, 
the correlation effect of crowded ions, the screening effect of polar water, as well as the 
charge/space competition effect of ions and water molecules of different sizes and valences. 
These effects are all closely related to the interstitial voids between particles and described 
by two additional terms, namely, the correlation length and the steric potential. 


4. A MOLECULAR-CONTINUUM MODEL OF A CALCIUM CHANNEL 


To test the PNPF theory, we use the Lipkind-Fozzard molecular model [3^ of L-type Ca 
channels in which the EEEE locus (four glutamate side chains modeled by 8 0^/^“ ions) 
forms a high-affinity Ca^’*' binding site that is essential to Ca^’*' selectivity, blockage, and 
permeation. We refer to Fig. 9 in or Fig. 1 in [l^ for a 3D illustration of the EEEE 
locus. A 2D cross section of a simplified 3D channel geometry for the present work is shown 
in Fig. 1, where the central circle denotes the binding site, the other four circles denote the 
side view of 8 0 ^^^“ ions, is the solvent region consisting of two baths and the channel 
pore including the binding domain UBind, dfls is the solvent boundary, and ^Death is the 
outside and inside bath boundary. Fig. 2 is a sketch of the binding site and 0^/^“ ions, 
where is the distance between the center of a binding Ca^"*" ion and the center Cj of any 
0^/^“, and A is any point on the surface of the site. In our model, the 8 0^/^“ ions are not 
contained in the solvent region Particle species are indexed by 1, 2, 3, and 4 for Na"*", 
Ca^"*", CD, and H 2 O with radii Oi = = 0.95, 02 = aQ^ 2 + = 0.99, 03 = = 1.81, and 

04 = an^n = 1.4 A, respectively. 

In [ 1 ^ , we proposed an algebraic model for calculating the electrical potential 0^ and the 

steric potential 5^ in hlBind by using Coulomb’s law with the atomic structure of binding 
ion and atoms in a channel protein as shown in Fig. 2, without solving the Poisson-Fermi 
equation (9) in hlBind- The binding potential 0^ is then used as a Dirichlet boundary condition 
in flBind for solving the PF equation in the solvent region between the bath and binding 
boundary, i.e., in r 2 s\r 2 Bind, to obtain the potential profile 0(r) that connects 0^ in UBind to 
the potential D (or 10 ,) on the inside (or outside) bath boundary. 


The filter domain defined in [l^ is simply taken to be the binding domain UBind D this 
paper. The volume of this domain is an unknown variable Vb that changes with different 
charges in the binding site. We do not define an ad hoc filter for which its volume is fixed 
to one value in one implementation and possibly to another value in other implementation. 
We show that the variable binding volume Vb plays an essential role in determining the steric 
potential in and around the binding site and consequently the hydrophobicity of the EEEE 
locus under different bath conditions. 
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FIG. 1: A simplified Ca channel geometry with baths, pore, and binding site. The channel is 
placed in a cubic box with the length of each side being 40 A. The solvent region consists of 
two baths and the channel pore with the boundary The binding site llBind is contained in 
but the 0^/^“ ions are not in Qg- The outside and inside bath boundary is denoted by SllBath- 


y 



FIG. 2: The binding distance between the center of the binding Ca^'*' ion and the center Cj of the 
jth 01 / 2 - 

ion is denoted by (Iq^ for j = 1, • • ■ ,8. A is any point on the surface of the binding ion. 


The algebraic model 


14 is defined in G 


Bind 


and consists 


of the following equations 


( 0\ = VbCf exp(-/3i0b + 

< Ol = VbCf exp{-f32(l)b +S^) , (18) 

(Ol = VbCf exp(S'fe"'') 



Vb - viO\ - V2O2 - v^O\ 

-Sirs- 


( 19 ) 
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( 20 ) 


where 0\, O^, and O 4 denote the occnpancy nnmbers of Na"*", Ca^’*', and H 2 O in Vb, respec¬ 


tively, (pb and are average electrical and steric potentials, and \cj — A\ is the distance 
between A and Cj in Fig. 2. 

In this mean held, we allow Oj and O 2 (and hence the total charge 0^ezj^g^+ + 02eZ(j^2+) 
to vary continnonsly snbject to the condition on their snm Oj -|- O 2 = 1 in the binding 
volnme Vb- Eqs. (18) and (19) nniqnely determine the fonr nnknowns Vb, O 4 , 0^ and 
with Oj and O 2 being given. Eq. (20) nniqnely determines the locations (cj) of 8 ions 

(and thns the binding distance (Iq^ or cIq^' in Fig. 2) once 0^, is obtained. Note that the 

binding distance changes continnonsly with varying O 4 and O 2 bnt 0 ^ 

remains hxed, where the binding ion O^Na + O^Ca is a linear combination of Na"*" and Ca^’*'. 
Therefore, 0^/^“ ions are movable — the protein is hexible in onr model — as their locations 


Cj changes with varying 0\ and O 2 iJ]. Note that we change the probability notation Pi 


m [ijj to the occnpancy notation Oi to rehect the deterministic, instead of probabilistic, 
natnre of the PNPF continnnm model. In this simple algebraic model, we do not consider 
the hydrogen ions that may react with carboxyl anions in the protein. Experiments done at 
pH 8 (as many have been done) do not involve association of hydrogen ions with carboxyl 
anions. 

For the half-blockage experimental condition 


= Cf = 32 mM, C^^ 2 + = Cf = 0.9 fiM, 

" -V-' 

Experimental Data 


( 21 ) 


we follow convention and assnme relative occnpancies of a hlled channel, 0\ = 0.5 and 
O 2 = 0.5, and thereby obtain 0^ = —10.48 ksT/e, = —1.83, and Vb = 4.56 14 . 

The binding experiments nsed a hxed (7® + = (Pf = 32 mM and varions Ca^’*' bath 


concentrations C^g 2 + = C'f that imply different 0\ and O 2 of Na^ and Ca^’*' occnpying the 
binding site. The occupancy numbers Oj and O 2 are determined by 


91 


^ qP = exp(-(ft - 


( 22 ) 


where 0^ was just obtained from the case of equal occupancy. The occupancy ratio in (22) 
thus deviates from unity as (7® is varied along the horizontal axis of the binding curve from 


its midpoint value Cf = 0.9 fiM as shown in Fig. 5 in |14 . 

Keeping (pb hxed is equivalent to assuming that the relation (22) between the occupancy 


and bath concentration ratios is linear [1^. Moreover, keeping 0^ hxed in (20) is equivalent 


to assuming that the 0 ^/^“ ions (cj) move continuously in response to the continuous change 
of charges 0'lz^g^+ -|- 02 ^ca 2 + iii the binding site. In [1^, the charge change from (Na"*" 
occupying the site) to Zq^ 2 + (Ca^+ occupying the site) rehects a change of pore radius of 
about 2.3 A that is surprisingly close to the value of 2 A obtained by MD simulations js^. 
Note that the vacuum permittivity eo is chosen in (20) since both MD models in 35|,|36| treat 

OV2- 

ions explicitly as shown in Fig. 1 (or Fig. 9 in (35|)- The Coulomb forces between 
the binding ion and 0 ^/^“ ions should therefore be calculated in vacuum since nothing is 
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in between these ions. Our numerical results can thus be verihed with those of MD. Of 
course, our assumptions in the linear model should be modihed if more accurate structural 
information can be used to provide an extra equation for variable 0^. The permittivity is 
chosen as 30eo [l^ in our forthcoming studies on the protein structure of a sodium/calcium 
exchanger (NCX) j^. Meanwhile, the linear model seems at least as good as the homology 
structure itself, and provides potentially useful and interesting insights as we show in the 
results section. Nevertheless, we imagine that nature might design flexible proteins so that 
is hxed or slightly perturbed by small thermal variations so that the linear model is still 
tolerable within numerical errors in theoretical simulations. 

The simple atomic structure in Fig. 2 elucidates algebraic and subsequent PNPF calcu¬ 
lations in a concise way. The molecular-continuum model presented here can be extended to 
deal with more complex nonequilibrium systems in real protein channels in future studies. 
Application of the algebraic model to the NCX structure is briefly discussed in 141. It 


will be interesting to apply the present model to recent structures of a TRPVl ion channel 
and a voltage gated calcium channel [^. 

For nonequilibrium cases, the binding steric potential Sf^ is assigned its equilibrium value 
in subsequent PNPF calculations, i.e., the void fraction F(r) in fleind is assumed to remain 
unchanged from equilibrium to nonequilibrium. The electrical potential 0^ will be modihed 


by the membrane potential Vi — Vo [70| and then used as a Dirichlet type condition for the 


potential function 0(r) in fleind- We summarize the boundary conditions for the PF (9) and 
NPF (15) equations dehned in the solvent region in Fig. 1 as 


0(r) = 0b(r) in flsind, 0(r) = K,i on Sflsath, 

V0(r) ■ n = 0 on d^ls\d^f, 

c,(r) = Cf = [Na+1„_„ C 2 (r) = Cj = |Ca"+]„,, Calr) = Cf = [Crl„_, on 
Ji(r) ■ n = 0 on 0fls\cXlBath, 


(23) 


where n is an outward unit vector on the solvent boundary dQg- Note that the electrostatic 
potential 0(r) is prescribed as a Dirichlet function 0fe(r) whose spatial average in Deind is the 
constant 0^. However, the binding domain flBind is treated as an interior domain instead of 
boundary domain for the hux equation (15). An iterative process of solving PF (9) and then 
NPF (15) is repeated until self-consistent solutions of 0(r) and C'i(r) are reached within a 
tolerable error bound. 

Treating the interior domain flBind C Ds in place of the conventional boundary d^ls for 
the potential condition 0(r) = 0fc(r) is not a conventional way to solve the Poisson equation. 
This method is needed because the binding potential 0^ determined by (18)-(20) is coupled 
to the steric potential S^, that in turn depends crucially on the conformation of the binding 
ion and protein atoms, their interstitial voids, and their charges as shown in Fig. 2. Water 
also plays a vital role as in (19). The steric equation (19) is in the interior domain flBind- 
In other words, for calculating S^, we need to consider voids and water volume that are 
interior quantities in flBind and cannot be specihed on the solvent boundary dfls- We thus 
need to impose 0(r) = 0fe(r) in flBind because 0fe(r) is given by 0^ in flBind- 

If a conventional method is used to solve the Poisson (or PF) equation [l3[, the resulting 
steric potential S'f, (as an output of 0(r) by (7)) may be completely incorrect in flBind 
because the atomic equations (19) and (20) are not used. We do not have any differential 
equation for the steric function 5*’''^(r) for which appropriate boundary conditions near flBind 
can be imposed if a conventional method is used. Moreover, the important role of water and 
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its volume effect is not taken into account in conventional methods or models. 

The models and methods proposed in this paper are still coarse approximations to ion 
transport as the PNPF theory is in its early development. Nevertheless, the theory proposed 
here provides many atomic properties such as (19) and (20) that have been shown to be 
important for studying the binding mechanism in Cay channels [l^ and are also important 
for the transport mechanism as shown in the next section. Incorporating atomic properties 
into continuum models is a step forward to improve and rehne the continuum theory that 
has been challenged for its accuracy when compared to (mostly equilibrium calculations) 
MC, BD, or MD j^, [s^, 71, 72\. Continuum models on the other hand have substantial 
advantages in efficiency that are of great importance in studying a range of conditions and 
concentrations, as are present in experiments. 


5. RESULTS 


Cay channel conducts primarily Na"*" when the Ca^’*' concentration is below 1 /iM but it 
conducts primarily Ca^"*" in the physiological concentration range mM. In j^, 19 extracellular 
solutions and 3 intracellular solutions were studied experimentally. The range of [Ca^’'']o is 
10®-fold from 10“^°'^ to 10“^ M as given in [^. Explaining the biological function of trace 
Ca^’*' concentrations is a crucial task of biophysical models while dealing with the large 
Ca^"*" concentrations found in extracellular solutions of all biological systems. This range 
of calcium concentrations poses severe obstacles for MD and BD simulations even on the 
most advanced computers to date J73( | . To our knowledge, a comparison of MD simulations 
with experimental measurements has not yet been reported without invoking arbitrary 
(i.e., untested) extrapolation methods for handling the 10®-fold variation of concentration 
and the dynamics of ionic flow BQ- 

PNPF results are in accord with the experimental data in as shown in Fig. 3 under only 
the same salt conditions of NaCl and CaCl 2 in pure water, without considering protons and 
other bulk salts in experimental solutions. With [Na’'']i = [Na’'']o = 32 mM and [Ca^’'']i = 0, 
the membrane potential is hxed at —20 mV (W = 0 and V = ~20 mV) throughout, as 
assumed in Fig. 11 of for all single-channel currents (in femto ampere fA) recorded in 
the experiment. Note that the experimental currents have been converted to single-channel 
currents in Fig. 11 of as shown on the right hand y-axis in that hgure. The currents 
are on the femtoscale because calcium channels have been long recognized to be in some 
sense blocked sodium channels and here too we have small membrane potentials. The small 
circles in Fig. 3 denote the estimated currents (by eye) from Fig. 11 of and the plus 
sign denotes the current calculated by PNPF. The channel is almost blocked by Ca^’*' ions 
with a current of about 15 fA at [Ca^’'']o = 10“^'^ M. Half-blockage current (about 77 fA) 
is dehned by the one half of the saturation current (about 154 fA at [Ca^’'']o = 10“^°'^ M). 
The half-blockage Ca^’*' concentration is about [Ca^’'“]o = 0.9 /iM and that is used to dehne 
the midpoint binding condition (21). 

PNPF deals naturally with the main experimental data of ionic flow based on this binding 
dehnition. Physical parameters in (9) and (16) are summarized in Table 1 with their physical 
meaning, numerical values, and units for ease of reference. The diffusion coefficient in the 
channel pore is taken as 6iDi for each ionic species, where Di are bulk values in the table 
and 9i = 0.1 are factors for the pore. All physical parameters are kept hxed throughout. 

PNPF can provide more physical details of ion transport inside the channel pore such as 
electrical potential (0(r) in Fig. 4), steric potential (S'*'’‘^(r) in Fig. 5), energy wells (£’ 2 ( 1 ') 







15 



-7-6-5 4 

log^gfCaCy in M added to 32 mM NaCf 


FIG. 3: Anomalous Mole Fraction Effect. Single-channel inward current in femto ampere (fA 
plotted as a function of logiQ[Ca^''']o. Experimental data marked by small circles are those in 
whereas the PNPF data are denoted by the plus sign. 


in Fig. 6), water density (C' 4 (r) shown in Fig. 7), void volume fraction (F(r) in Fig. 8), 
dielectric function (F(r) in Fig. 9), crowded ions in binding site (Fig. 10), concentrations 
(G 2 (r) in Fig. 11), and flux densities (|J 2 (r)| in Fig. 12). The electric potential prohles 
remain almost unchanged for various (7® as shown in Fig. 4 following the linear model of 
the occupancy equation (22). 

What is new? The steric potential prohles shown in Fig. 5 represent the novelty of 
the PNPF theory. All effects of volume exclusion, interstitial void, conhguration entropy, 
short range interactions, correlation, polarization, screening, and dielectric response of this 
nonideal system are described by the steric functional S''^'^'^(r). The Ca^^ energy landscape 
^^ 2 ( 1 ') = (/920(r) — S^'^^{r))kBT (Fig. 6) is modihed by this potential. The steric potential 
in the binding region decreases drastically from —1.30 to —10.34 ksT as increases from 
10“’^-^ to 10“^ M. The water density C' 4 (r) = Gf exp(S'*''‘^(r)) (Fig. 7) and the void volume 
F(r) (Fig. 8) decrease as well because more Ca^’*' ions in the bath make the binding site 
more crowded as was previously seen in the algebraic Fermi model fl3 |. 
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FIG. 4: The electrical potential (/)(r) profiles (averaged over each cross section) along the pore axis 
for various Cf ranging from M to 10“^ M. All the following figures are obtained with the 

same averaging method and the same range of Cf. 


Table 1. Notations and Physical Constants 


Symbol 

Meaning 

Valne 

Unit 

ks 

Boltzmann constant 

1.38 X 10-23 

J/K 

T 

temperatnre 

298.15 

K 

e 

proton charge 

1.602 X 10-19 

C 

eo 

permittivity of vacnnm 

8.85 X 10-1^ 

F/cm 

6w 

water dielectric constant 

78.5 


Ic 

correlation length 

1.98 

A 

Di 

Na"*" diffnsion coefficient 

1.334 X 10-3 

cm2/s 

D2 

Ca^"*" diffnsion coefficient 

0.792 X 10-3 

cm2/s 

Ds 

Cl~ diffnsion coefficient 

2.032 X 10-3 

cm2/s 

Cf 

Na"*" bath concentration 

32 

mM 


Ca^’*' bath concentration 

10-10-3 ~ 10-2 

M 


inside (ontside) voltage 

0 (-20) 

mV 


In physiological bath conditions = 10“^ ~ 10“^ M, Fig. 7 shows that the region 
containing the binding site with the length abont 10 A is very dry (hydrophobic), which 
agrees with the recent crystallographic analysis of the Ca^’*' binding site of the related 
protein NCX with the EETT locns showing a hydrophobic patch (also abont 10 A in length) 
formed by the conserved Pro residnes. The hydrophobicity near the binding site in onr 
model is described by the continnons water density fnnction C' 4 (r) via the continnons steric 
fnnction S'*’''^(r) as shown in Fig. 5. At Cf = 10“^ M, the magnitnde of the steric energy 
= —10.34 ksT is comparable to that of the electrostatic energy cj) = —10.48 ksT/e. This 
snrprisingly large energy dne only to the steric effect has not been qnantihed and observed 
by MD, MC, or other continnnm methods in Cay channel modeling, as far as we know. As 
observed from Fig. 3, ionic transport is blocked by the competition between Na"*" and Ca^"*" 
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FIG. 5: The averaged steric potential profiles. 



FIG. 6: The averaged Ca^^ energy well £^ 2 (r) profiles. 


ions in the range = 10“^'^ ~ 10“^'^ M. In this blocking range, the corresponding steric 
prohles in Fig. 5 are wider indicating that the water density or the void volnme is more 
evenly distributed. 

Fig. 9 demonstrates the combined effects of correlation, polarization, and screening in this 
highly inhomogeneous electrolyte by means of the variation of dielectric coefficient produced 
by the Poisson-Fermi equation (9). Note that the dielectric function can also be calculated 
by ?(r) = Cb + G 4 (r)(ew — €-b)/C^ using the water density function G 4 (r) as proposed in j75j |. 
where = 2. Fig. 10 illustrates how ions crowd in the highly charged binding site and, in 
the meantime, voids and water vacate under the condition G® = 10“®'^ M. 

The Ca^’*' occupancy 0\ in the binding site increases from 0.69 (not shown) at G® = 10“®'^ 
M to almost 1 at G® = 10“^ M. The corresponding peak Ca^’*' concentration shown in Fig. 
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FIG. 7: The averaged water density C 4 (r) profiles. 



FIG. 8: The averaged void volume fraction F(r) profiles. 


11 also increases from 261.83 to 408.53 M in this range. The largest concentration is still 
below the maximum allowable value = l/v 2 = 408.57 M as implied by the Fermi 

distribution (7). We obtained this incredibly large concentration since 0\ k. 1 vn (19), 
which in turn yields l/ub ~ 408.53 M as ~ V 2 indicating the importance of the atomic 
nature of the steric potential (19) with variable v^. This hgure demonstrates that the PNPF 
model can capture the atomic properties of ions in flow, a critical (however small) step in 
continuum theory toward the ultimate accuracy in theoretical simulations. From Fig. 3, we 
observe that Ca^^ currents increase dramatically when [Ca^’'']o increases from 10 “^'^ to 10 “^ 
M in the physiological mM range of Cay channels. The corresponding flux density prohles 
in the binding region increase dramatically too as shown in Fig. 12. 

We make a hnal remark about these results. As intensively studied in 113, HQ , the 
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FIG. 9: The averaged dielectric function e(r) profiles. 



FIG. 10: Crowded Charge in Binding Site. The volume fractions of voids, water, and ions per unit 
volume at the Ca^'*' bath concentration C® = 10“®''^ M. 


classical PNP model fails to yield ionic concentrations and currents properly (compared with 
those of Brownian dynamics), especially for narrow channels, because the classical model 
does not include the volume effect of ions, water, and voids, the correlation effect of ions, and 
the screening effect of water. The PNPF model not only computes ionic currents comparable 
to the experimental data (Fig. 3) for the narrow calcium channel but also provides many 
physical properties (Figs. 5-10, for example) that are shown to depend critically on the 
proposed steric potential. The PNPF model overcomes the limitation of the classical PNP 
model with respect to these effects and properties. 
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FIG. 11: The averaged Ca^'*' concentration C' 2 (r) profiles. 



FIG. 12: The averaged Ca^"'' flux density |J 2 (r)| profiles. 


6. CONCLUSION 

We propose a Poisson-Nernst-Planck-Fermi model for studying equilibrium and nonequi¬ 
librium systems of ionic liquids and solution electrolytes. The excluded volume effect of 
different sizes of ions and water molecules, the correlation effect of crowded ions, and the 
screening effect of polar water molecules in inhomogeneous aqueous electrolytes are all in¬ 
cluded in this model. The model was verihed by a set of experimental currents of L-type 
Ca channels recorded in a 10®-fold range of Ca^"*" concentrations that show the exceptional 
selectivity of these channels. 

We also propose a consistent Gibbs free energy functional leading to a Fermi like dis¬ 
tribution of hard spherical particles in the electrolytic system. The Gibbs-Fermi functional 
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is shown to converge to the Gibbs-Boltzmann functional that yields the Boltzmann distri¬ 
bution as the volumes of all particles and the correlation length approach zero. Moreover, 
we introduce a Gibbs-Fermi entropy for which the excluded volume of water molecules and 
the dynamic distribution of interstitial voids between particles are needed to establish a 
consistent generalization of nonideal inhomogeneous electrolytes for both equilibrium and 
nonequilibrium systems. The local probability of any particle species in the Gibbs-Fermi 
entropy is expressed in terms of the local and bulk ratios between the particle and void frac¬ 
tions per unit volume. We show that if all particles are treated as hard spheres, the voids 
must be included in the Gibbs free energy functional. The voids are created by packing 
defects and Lennard-Jones and Goulomb forces between particles. The void effect plays an 
essential role in our theory as well as in our results, especially in the steric energy. 

Most of the results in this article seem to be novel, because consistent models including 
voids, water volume, and Fermi distributions have not been developed previously, as far 
as we know. These numerical results provide useful tools to develop insight into a vari¬ 
ety of physical mechanisms ranging from binding, to permeation, blocking, flexibility, and 
charge/space competition of the channel. 
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